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NONPARAMETRIC  ESTIMATION  OF  THE  PROBABILITY  OF 
A  LONG  DELAY  IN  THE  M/G/1  QUEUE 

BY 

D.  P.  GAVER 
P.  A.  JACOBS 

OPERATIONS  RESEARCH  DEPARTMENT 
•  NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA   939^3-5000 

1  .   Introduction 

The  application  of  probability  theory  to  a  wide  variety  of  congestion 
problems  has  been  well  catalogued  in  many  papers.  Most  of  the  elegant 
solutions  obtained  are  in  somewhat  implicit  form,  being  presented  as 
functional  equations,  or,  frequently,  as  integral  (Laplace)  transforms, 
generating  functions,  and  sometimes  as  combinations  of  the  above.  The 
results  obtained  naturally  appear  in  terms  of  component  distribution 
functions  and  stochastic  processes  (renewal,  Poisson  etc).   Only  rarely  are 
issues  addressed  that  arise  when  actual  data  is  to  be  used  as  a  basis  for 
inference  from  the  models;  however,  see  Cox  [1965]. 

In  this  paper  we  consider  the  nonparametric  estimation  of  the 
probability  of  a  long  customer  delay  in  an  M/G/1  system,  given  a  known 
Poisson  arrival  rate  A  and  observations  of  independent  service  times  from 
the  service  distribution,  presumed  unknown.   Although  the  approach  and 
results  are  given  concretely  for  the  M/G/1  system,  they  apply  more  widely. 

To  be  specific,  consider  a  single  server  system  approached  by  a 
stationary  Poisson  (A)  traffic  with  A  known.   Service  times,  X  are 
independent  identically  distributed  with  unknown  distribution  function  F  ; 
assume  AE[X]  =  p  <  1 .   Let  observations  of  the  service  times  be  all  that  is 


known  about  F;  denote  these  by  x  ,  x  ,  ...,  x  .   The  objective  is  to  supply 
nonparametric  estimates,  and  error  assessments  thereof,  of  the  probability 
of  a  long  delay  experienced  by  an  arriving  customer. 

It  is  well  known  that  if  W(t)  is  the  virtual  waiting  time  in  the  M/G/1 
queue  and  p  <  1 ,  then  the  moment  generating  function 


nr  SWi    ,.   nr  SW(t)i  , ,  , x 

E[e   J  =  lira  E[e     J  (1.1) 


=  (1-p)[1  -  pA(s)]"'1 


where  A(s)  is  the  moment  generating  function  of  a  distribution  H.   If  A(s) 
exists  for  s  <  s  ,  s  >  0,  then  there  will  be  a  smallest  real  zero  s  =  <  >  0 
of  the  denominator  of  (1.1)  which  can  be  used  to  show  that 

P{W  >  w}  -  D(<)  e"KW,        w  -  «.  (1.2) 


We  will  always  assume  this  is  the  case. 

One  way  of  establishing  (1.2)  is  to  introduce 

sW         °° 
4,(3)  =  E[e ^_1I  =  f  p{W  >  w}  e3Wdw  (1.3) 


=  f  P{W  >  w 


into  (1.1)  and  to  rewrite  in  the  form 


4,(S)  =  p[A(3)  "  1]  +  PA(s)f(s)  (1.4 


-2- 


which  is  equivalent  to  the  terminating  renewal  equation,  see  Feller  [1966] 
p.  362, 


Fw(w)  =  P{W  >  w}  =  PH(w)  +  p  f  P{w  >  w-x}  H(dx)  (1.5) 

0 


where 

0 


-  #      -      0w 
Following  Feller,  introduce  F   (w)  =  F  (w)  e   ,  9  real  and  positive,  into 

w       w 

(1.5)  to  obtain 


w 
Fw#(w)  =  p  H#(w)  +   J  F#(w-x)  pe9XH(dx).  (1.7) 

0 


Choose  0  =  <  so  that 


i  H(dx)  = 

0 


[  pe9xH(dx)  =  [  H#(dx)  =  1  (1.8) 


yielding  a  standard  renewal  equation  for  F  .   From  the  key  renewal  theorem, 

w 

it   follows   that 


lim  F,,(w)    =   lim  F,.(w)e        =      ,    >.  (1.9) 

W  W  m( < ) 

W->°°  W  +  00 
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where 


c(<)  =  P  f  eKX  H(x)dx  (1.10) 


and 


m(<)  =  p  f  xe   H(dx).  (1.11) 

0 

Summarizing 


P{W  >  t}  -  ^^■   e  Kt   t  ■*   «,  (1.12) 


m(<) 


where  <  is  the  positive  solution  to  the  equation 


X9  1U(0)  -  1]  =  1  (1.13) 


with 


.(9)  =  E[e9X].  (1.U) 


In  section  2,  a  nonparametric  estimate,  <,  of  <  is  studied  which  solves 
equation  (1.13)  with  the  empirical  moment  generating  function  replacing  <j>. 
This  estimate  is  related  to  an  estimate  studied  by  Stigler  [1971]  in  the 
context  of  estimating  the  probability  of  extinction  for  branching  processes. 
k   is  shown  to  be  a  consistent  estimate  of  <  having  a  distribution  in  the 
domain  of  attraction  of  a  stable  law.   Under  certain  conditions,  a  central 
limit  theorem  for  <   can  be  obtained. 
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In  section  3  a  non parametric  estimate  of 


P(t)  s  ^T  e_Kt  (,-15> 


is  given. 

In  Section  4,  some  results  of  simulation  studies  of  the  estimates  of  < 
and  p(t)  are  presented. 

2.   Nonparametric  Estimation  of  the  Exponent  <  of  the  Probability  of  a 

Long  Wait 

Assume  for   the  remainder   of   the   paper   that    A=1    is   known.      Let 

x„  ,    ....    x     be   the  observations   of   the  service   times.      A  non-parametric 
1  n 

estimate  of   the  moment   generating  function  of   the   service   time  distribution 
is 


n        9x. 
.(9)    =  -     I     e      l.  (2.1) 

n   i-1 


The  sample  equivalent  to  equation  (1.13)  is  thus 


1  =  e  1[J(6)-1].  (2.2: 


At  e  =  0,  the  RHS  of  (2.2)  is  x  which  is  less  than  1  if 

X   +  X   +  ...  +  X 

x  =  — <  1  ;  the  data  can  only  be  analyzed  for  a  stationary 

n 

model  under  this  assumption,  which  will  be  made  in  what  follows.   Further, 
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(2.2)  has  a  unique  solution,  k,  which  is  an  estimate  of  k.   A  three-term 
Taylor's  expansion  of  the  RHS  of  (2.2)  about  9=0  gives  an  initial  estimate 


<H  =  2[1  -  x]  (x2)"1  (2.3) 


where  x  =  —  [x,  +  ...  +  x  ],  the  k  sample  moment.   Since 
n   1         n 

i  *  —   i   __  "  ■*■ 

e"1  [<J)(9)  H ■  1  ]  ^  x  +  t;   9  x2,  <  is  an  upper  bound  for  <.   Equation  (2.2)  can 

d  H 

>\ 

be  solved  via  search  or  Newton-Raphson  iteration  to  obtain  the  estimate  k. 
We  will  now  present  asymptotic  results  for  the  distribution  of  <  as  the 

sample  size  n  ■*  ».   By  assumption,  if  X  is  a  random  variable  having  the 

B'X 
service  time  distribution,  then  E[e   ]  <  *  for  some  B'  >  <  >  0.   Thus, 

there  is  B  >  <  such  that  for  all  0  <  b  ^  B 


E[XnebX]  <  oo.  (2. it) 


It  follows  from  the  monotonicity  of  the  RHS  of  (2.2)  that 


PU  >  B)  -  PJttBJ  «  1  „  [»(B)  -  1]  <    ,  .  U(B)  -  1]  , 

B  D  D 

=  Pji(B)  -  <|,(B)  <  B[1  -    U(Bg  -  1]]|.  (2.5) 


Since  4>(0)  is  a  monotone  function  and  B  >  <,  [1  -  B  [<J>(B)  -1]]  is  negative 
Thus,  by  the  strong  law  of  large  numbers 
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lim  P{<   >  B}  =  0.  (2.6) 

n-K*> 


Let 


f(e)  =  1  -  e  1[<t>(e)  -  1].  (2.7) 


Expand  f(<)  in  Taylor's  series  about  the  solution  <  of  (1.13).   Since 
f(<)  =  0, 

0  =  f(<)  +  (<-<)  f(<)  +  ^  (<-<)2  f"(6<)  (2.8) 

for  some  e.  Thus 


k  -  <  =  -f(<)[f'(<)  +  j   (<-<)  f"(6<)]  1;  (2.9: 


E[f(<)]  =  0; 


1  k'Y      k'Y 

E[f'(<)]  =  -^  E[-<Xe   +  eK     -  1]  =  3  <  0.  (2.10) 


By  the  strong  law  of  large  numbers 


lim  f(<)  =  0  (2.11) 

n-*-°° 

and 

lim  f'(<)  =  6  <  0  (2.12) 

n->-«> 
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with  probability  1.   It  follows  from  (2.6),  (2.9),  (2.11)  and  (2.12)  that  < 
converges  in  probability  to  <   as  n  ■*  °°.   Thus,  <  is  a  consistent  estimate  of 


KC. 


If  E[e   ]  <  «,  then 


Var[fU)]    -  -  -I  Var[eKX]    =  -  r2.  (2.13) 

n   <■*  n 


By  the   central   limit   theorem 


£  /~n~   {<  -  <}    -  £  /~n~  [-f(<)]    [f'(tc)    +  1  (<-<)    f"(9<)]"1  (2.14) 


is  asymptotically  standard  normal  as  n  -*■   °°. 

2<X  <X 

If  E[e   ]  =  °°,  then  the  distribution  G  of  e  L  is  in  the  domain  of 

attraction  of  a  stable  law  with  index  1  <  a  <  2,  see  Feller  [1966],   Let 

{a  }  be  a  sequence  of  numbers  such  that 
n 


a 
\  [n  x2  G(dx)  +   1        as  n  ->  «.  (2.15) 

*n  0 


The  normalized  random  variable 


n   kX 

—  [<-<]=  —  [--  I     e   i][f'(<)  +  U<  -   <)  f"(6<)]_1       (2.16) 
a  an..  2 

n  n    i  =  1 


is  in  the  domain  of  a  stable  law  with  index  a,  where  1  <  a  <  2  is  such  that 

x  G(dx)   -  y     L(y)  where  L  is  a  slowly  varying  function. 
0 

To  summarize,  we  have  the  following  result. 


PROPOSITION        a.   <  is  a  consistent  estimate  of  k. 

2kX 

b.  If  E[e   ]  <  °°,  then  the  distribution  of  <-<  is 

asymptotically  normal. 

2<X 

c.  If  E[e   ]  =  °°,  then  <-<  is  in  the  domain  of 

attraction  of  a  stable  law. 


We  will  suppose  for  the  remainder  of  this  section  that  the  service  time 
distribution  is  exponential  with  mean  —  >  — .   In  this  case,  (since  A  =  1  by 
assumption) , 


<  =  u  -  1 ;  (2.17) 

Var[e<X]  =  m(m-1)2  (2-y)"1 ;  (2.18) 

and 

E[f(<)]  -  -1.  (2.19) 

From  the  central  limit  theorem  (2.14) 


nVar(<)  -  (§]2  =  u(2-u)  1  =  f^l.  (2.20) 


-9- 


2   1/2  -1 

Therefore,    if   g(x)    =   (1-x    )  -    1    -   sin      (-x),    then 


nVar[g(<)]    -   [g'(<)]   nVar[<]    -   1 


Thus,  g  is  an  asymptotically  variance-stabilizing  transform  of  <.  The 
simpler  related  transformation  ln(1+<)  is  used  in  the  simulations  of  < 
reported   in  Section   4.    ' 

3.      A  Nonparametric  Estimate  of   P{W  >   t} 
The   asymptotic  analytic  result   is 


P{W  >   t}    -  ||^  e  Kt  =  p(t)  (3.1) 


as   t  ■*  °° 


where 


c(<)    =    [  e<y      f  P{X  >   x}dx  dy  (3-2) 


2 


=   \     Y~  PUedx}    +  ^t     |   Ee<X  -   1   -   kx  -  -^-]    P{Xedx} 


and 


m(<)    =      I"   xe<X  P{X  >   x}   dx  (3-3) 
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1 

— r  1  I  Kxe 


J  <xe<X  P{Xedx}  -  f  (e<X  -  1)  P{Xedx}} 


0 


=  1  |   x2   P{Xedx}    +  1  <    f   x3   P{X£dx: 
0  0 

CO 

+    I"    (kx  -    1  )    R(x)    P{Xedx} 
0 


where 


(<x)' 


n/    v         r    <x  .                           ^<x;    l 
R(x)    =    [e        -    1    -    kx    -   r J 


(3.4) 


An   estimate  of   P{W  >   t}    is 


*/«.  \        c(<)         -<t 
p(t )    =  ^— ^ —  e 

m(<) 


(3.5) 


where  <  is  the  positive  solution  to  equation  (2.2); 


c(<)  =  1  71  +  i-  R1(<); 
<2 


(3.6) 


m(<)  =  -  x2  +  p  <  x3  +  R2(<): 


(3.7) 


n    kx  . 


(kx.  )' 


R/k)  -  1  I     (e  X  -  1  -  <x<    -ji-  ); 
i-1 


1     2 


(3.8) 
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2 
~  ^     ^    n      <x .       ^     (kx.) 

R2(<)  -  (-  R1(<))  +  k  -     I     x.(e  l  -  1  -  kx. j—   ).      (3.9) 

i  =  1 


The  forms  of  c(<)  and  m(<)  are  chosen  for  the  numerical  stability  of  the 
ratio  c(<)  m(<) 

In  the  remainder  of  this  section  we  will  assume  the  service  time 

1  1 

distribution  is  exponential  with  mean  —  >  —  and  x  <  1.   In  this  M/M/1  queue 

M  2 

case,  it  is  well  known  that 


p(t)  =  -  exp{-(M-1)t}  =  (1  +  <)  1  e  <t  (3.10) 


from  (2.16). 

We  will  now  motivate  a  transformation  of  p(t)  which  is  used  in  the 

y 
simulation  studies.  Let  Y  =  ln(1+<),  then  <  =  e  -  1  and 


p(t)  =  exp{-Y  -  [eY  -  1]t}.  (3.1D 


Let 


h(Y)  =  In  p(t)  =  -  Y  -  [eY  -  1]t  (3.12) 


and   Y   =   ln(1    +   k)  .      A  Taylor's   expansion  yields 


nVar[Y]   =  nVar[ln(1    +   <)]  (3-13) 


_1_ 
(1+k) 


=  n   , ,  ,    w  Var[<] 


1-K 


1        ■   [1-(eY-   I)2]"' 


-^2 
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from   (2.20).      Hence, 


nVar[h(Y)]    =    [-1    -e\]2    [1    -    (eY  -    I)2]    1 

=    [-    1    -   t   -    (eY  -    1)t]2[1    -    (eY  -    I)2]"1.  (3.1U) 


It  follows  from  the  definition  of  h(Y)  that 


nVar[h(Y)]  =  [-(1+t)  +  h(Y)  +  Y]2{1  -  (eY  -  1)2}  1  .      (3-15: 


Since  0  <  k  =  eY  -  1  <  1  ,   0  <  T  <  In  2.   Thus,  for  t  large 


nVarCh(Y)]  =  [-h(Y)  +  t]2  (3-16) 


which  suggests  that  the  transformation  ln[t  -  In  p(t)]  will  tend  to 
stabilize  the  variance  of  p(t),  at  least  for  exponential  service  time 
distributions. 

4.   Simulation  Results 

In  this  section  results  are  presented  of  a  simulation  experiment  to 
study  small  sample  behavior  of  the  estimates  of  <  and  P{W  >  t}.   All 
simulations  were  carried  out  on  an  IBM  3033AP  computer  at  the  Naval 
Postgraduate  School  using  the  LLRANDOM  II  random  number  generating  package; 
Lewis  and  Uribe  [1981].   In  each  replication  50  exponential  service  times 
are  generated.   If  the  mean  of  the  service  times  is  larger  than  1,  equation 
(2.2)  has  no  positive  solution.   In  this  case,  another  50  service  times  are 
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generated.      The  estimates,    <  of    (2.2)    and   p(t)   of   equation    (3*5)    are 
computed  for  each  replication. 

Standard  errors   of   the   estimates   are  estimated   in  two  ways.      One   is   the 
appropriate   asymptotic   variance  expression.      The  other   is   the  jackknife 
procedure. 

The  jackknife   is  a  procedure  originally   introduced  by  Quenouille  for 
bias  reduction   [1956],    and  adapted   by  Tukey   [1958]   to  obtain  approximate 
confidence   intervals.      Suppose   interest   is   on  a  parameter    8    (e.g.    k  or   p(t)) 
that    is  estimated   by    0  using  a   complex  calculation  from  data  x    ,...,x    .      The 
idea  is  that  of  assessing  variability  by  recomputing  0  after  removing 
independent   subgroups   of   data  of   equal   size,   and  then   using  the  recomputed   0 
values  to  estimate  a  variance  which  is  in  turn  applied  to  state  a  standard 
error   or   a  two-sided  confidence   interval   that   contains   the  true   0  with 
specified  confidence.      A  few  more  details  follow;    for  more,   see  Efron  [1982] 
and  his  more  recent   work,   or  Mosteller   and  Tukey   [1977].      The   actual 
calculation  involves   splitting  the  n   data   points   into   g  disjoint   groups   of 
size  m;    n=mg.      In  our   simulations   g  is  always    10.      Then   calculate   0/    .%» 
j=1,2,...,g:      the  estimate  of   0  based  on  a  reduced  data  set   that   omits   the 
j        group.      In  the  simulations,    the   first   group   is   the   first  five 
(unordered)   service  times;    the  second  group  is  the  second  five  service 
times,    etc. 

Now  Tukey  computes   pseudo-values 


yj   -go  -  (g-1)   0{.J5 
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which  are  treated  as  independent.   Tukey  recommends  referring  the  mean  of 

the  pseudo-values  y  to  Student's  t  with  g-1  degrees  of  freedom  to  obtain 

confidence  limits. 

In  the  jackknife  procedure  to  estimate  k,  it  is  sometimes  the  case  that 

a  positive  solution  of  equation  (2.2)  does  not  exist  for  the  data  set  that 

omits  the  j    subgroup  because  the  mean  of  this  reduced  data  set  is  larger 

than  1.   In  this  case,  <,    ..    is  set  equal  to  the  smallest  of  the  <,.  ,'s  that 

(-j)  (1) 

could  be  computed  for  the  sample  considered.. 

Similiarly,  in  a  jackknife  procedure  to  estimate  p(t),  it  is  sometimes 

the  case  that  either  <(_-\    does  not  exist  or  P/.s(t)  exceeds  1  for  the 

reduced  data  set  that  omits  the  1   subgroup.   In  this  case  p,  ,.  is  set 

(-J) 

equal  to  the  largest  of  those  p,_..'s  that  could  be  computed  and  were  less 
than  1  . 


4. 1   Results  of  a  Simulation  Experiment  to  Estimate  < 

In  this  subsection  results  are  given  of  a  simulation  experiment  to 
estimate  k.   For  each  replication,  three  different  estimates  of  <  were 
computed.   Estimate  I,  <  ,  is  computed  by  numerically  solving  (2.2)  by 
search  starting  with  the  initial  value  <u  of  (2.3).   Estimate  II,  < 

H  11 

is  obtained  by  jackknifing  <  using  ten  subgroups;  the  mean  of  the  pseudo- 
values  is  the  estimate.   Estimate  III  is  obtained  by  jackknifing  ln(l+<) 
using  ten  subgroups;  the  inverse  transform  of  the  mean  of  the  pseudo-values 

ey  -  1  is  used  as  the  estimate  of  k.   In  the  cases  of  Estimates  II  and  III, 
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if  the  estimate  is  negative,  it  is  set  equal  to  0.   For  each  estimate,  the 
average  bias 


1   500   . 

B  "  500  l    (vK) 

i  =  1 


and  the  relative  mean  square  error 


500  k .-<  2 
Hel  MSE  -  goo   I  I—] 

1=1 


are  computed  where 

k  =  y-1 
in  the  case  of  exponential  service  times  with  mean  —  <  1. 

Results  of  the  simulation  appear  in  Table  1 .   All  of  the  estimates  have 
about  the  same  relative  mean  square  error.  Jackknife  estimates  II  and  III 

have  smaller  bias  than  the  straightforward  estimate  I.   Jackknifing  <  itself 

1 
rather  than  ln(<+1)  gives  the  smallest  bias.   As  —  increases  all  the 

u 

estimates  have  increased  relative  mean  square  error. 

A  simulation  study  was  conducted  to  compare  the  performance  of 
different  confidence  interval  procedures.   For  each  replication  three  80^ 
confidence  intervals  were  constructed.   Interval  procedure  I  is  a  normal 
confidence  interval  which  uses  the  straightforward  estimate  of  <,  <  ,  as  the 
point  estimate;  the  estimate  of  the  variance  is  the  data  version  of  the 
asymptotic  variance  in  the  central  limit  theorem  (2.1*0; 
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~2  *2>2      1 


<x 


CLT  H        n-1 


I   (e      1   -   $(<))  (4.1) 

i 


where 


<x. 
>(<)    =  -  I  e      1  (4.2) 


n 

l 


6    =  -  I  e      l   -    1    -   <x.e      *  (4.3) 

i 


with  <  =   k      in  this  case.      The   805J-point   of   the  normal   distribution   is  used 
to  construct   the   interval.      Limits    that   are  negative   are  set   equal    to  0. 

Confidence  interval   procedure  II    is   a  jackknife   confidence   interval 
which  jackknifes   k  and   uses    the   80%  point   of   the   student   t-distribution  with 
9   degrees   of   freedom.      Limits   that   are  negative   are  set  equal   to   0. 

Confidence   interval    procedure   III    is  a  jackknife  confidence   interval 
which  jackknifes   ln(<+1 )    and   uses    the   80%   point   of   the  student 
t-distribution  with   9  degrees   of  freedom  to  give  a   confidence   interval   for 
ln(<+1).      The   inverse  transformation  of   the  endpoints   of   the   interval   gives 
an   interval    for   <;    limits   that   are  negative   are  set   equal    to  0. 

Results   of   the   confidence   interval    simulation  appear   in  Table   2. 
Reported  are  the  number   of   500    intervals   that   cover   the   true   value  of 
k,    (C);    the   number   of   the   500   intervals   such   that   the  entire   interval   lies 
below   k,    (B);    and  the  number   of   the   500   intervals   such   that   the  entire 
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interval  lies  above  the  true  value,  (H).   The  average  length  of  the 
confidence  intervals  is  also  given. 

The  number  of  intervals  that  cover  the  true  value  of  k  for  procedure 


III  is  within  .80  ±  (1.96)/  J^  (.2)(.8)  =  [.765,  .835].   All  but  one  case 

DUU 

of  the  normal  confidence  intervals  of  procedure  I  are  outside  this  range. 
All  but  one  case  using  confidence  interval  procedure  II  are  inside  this 
range.   The  average  width  of  confidence  intervals  for  procedures  II  and  III 
are  about  the  same.   Thus,  although  the  jackknife  estimate  <TTT  is  a  little 
more  biased  than  <  ,  the  coverage  of  the  jackknife  confidence  interval  for 
<    appears  to  be  somewhat  better  than  for  jackknife  confidence  interval 
procedure  II. 
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TABLE  1 
Bias  and  Mean  Relative  Square  Error 
for  Estimates  of  k 


Distribution 

Expor 

lential 

Exponential 

Exponential 

Exponential 

1 

0.6 

1-0.7 

-  =  0.8 
u 

1=0.9 

Estimate 

<  = 

.6667 

<   -  .4286 

<  =  .2500 

<  =  .1111 

B   Rel  MSE 

B   Rel  MSE 

B    Rel  MSE 

B   Rel  MSE 

(S.E) 

(S.E.) 

(S.E)   (S.E.) 

(S.E)   (S.E.) 

(S.E)   (S.E.) 

(!) 

(!) 

(!) 

(!) 

I 

.1147 

.2350 

.0557   .3292 

.0057   .6026 

.0865  2.400 

(.OH) 

(.021) 

(.011)  (.029) 

(.008)  (.054) 

(.0067)  (.2134) 

[.172] 

[.130] 

[.023] 

[.779] 

II 

.0219 

.2256 

.0067   .3095 

.0102   .5382 

.0469   1.785 

(.014) 

(.019) 

(.011)  (.025) 

(.008)  (.046) 

(.0063)   (.161) 

[.033] 

[.016] 

[.041] 

[.422] 

III 

.0533 

.2212 

.0160   .3073 

.0268   .5576 

.0608  1.997 

(.014) 

(.019) 

(.011)  (.027) 

(.008)  (.049) 

(.008)   (.177) 

[.080] 

[.037] 

[.1072] 

[.5472] 
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Table  2 
Confidence  Intervals  for  < 


Procedure 
Distribution 

Coverage 

Width 

I 
B    (C)    •  H 
%        (%)           % 

II 
B     (C)     H 
%           (%)           % 

III 
B     (C)    H 
%           (%)         % 

I      II     I] 

AV     AV     At 

(S.E.)   (S.E.)  (S. 

Exponential 

-  =  0.6 

k   =  .6667 

Exponential 
1-0.7 
k  =  .4286 

Exponential 

-  =  0.8 
M 

k   =  .2500 

Exponential 

1-0.9 
M 

K    =  .1111 

45  (363)   92 
9  (72.6)   18.4 

49   (380)   71 
9.8  (76)   14.2 

28   (407)   65 
5.6  (81 .4)  13 

1   (429)   70 
.002(85.8)  14 

38   (408)  54 
7.6   (81.6)  10.8 

50   (409)   41 
10   (81.8)   8.2 

38   (420)   42 
7.6   (84)   8.4 

54    (408)   38 
10.8   (81 .6)  7.6 

29   (395)  76 
5.8   (79)  15.2 

42   (399)   59 
8.4  (79.8)  11.8 

33   (413)  54 
6.6  (82.6)  10.8 

52   (392)  56 
10.4  (78.4)11.2 

.6402   .8089   .77 
(.005)  (.0124)  (.01 

.5279   .6255   .61 
(.004)   (.009)  (.OC 

.4441   .4831    .45 
(.005)  (.008)   (.0C 

.3733  -3763   .39 
(.005)   (.009)   (.0 

4.2  Simulation  Results  for  Estimating  p(t) 

In  this  subsection,  results  are  given  of  a  simulation  experiment  to 
estimate  p(t).   For  each  replication  three  estimates  of  p(t)  are  computed 
Estimate  I,  pT(t)  is  computed  for  each  replication  using  formulas  (2.2), 
(3.5)  -  (3.9).   If  PT(t)  exceeds  1,  it  is  set  equal  to  1. 
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There  are  (at  least  two)  possible  ways  to  implement  a  jackknife 

procedure  to  estimate  p(t).   Estimate  II  is  obtained  by  jackknifing  ln(<+1); 

an  estimate  of  <  is  obtained  by  the  inverse  transformation  of  the  mean  of 
the  pseudo-values  y   > 


yLK    , 
<u    -  e    -  1; 


<   is  used  in  formulas  (3.5)  -  (3-9)  to  obtain  the  estimate  p   (t).   If 
p   (t)  exceeds  1,  it  is  set  equal  to  1. 

Estimate  III  is  obtained  by  jackknifing  ln[t  -  In  p(t)];  if  p(t) 
exceeds  1  for  a  reduced  data  set  that  omits  the  j    subgroup,  the  estimate 
of  ln[t  -  In  p(t)]  for  that  reduced  data  set  is  put  equal  to  the  smallest 
estimate  that  could  be  computed  from  the  other  reduced  data  sets.   An 
estimate  of  p(t)  is  obtained  by  inverse  transformation  of  the  average  of  the 
pseudo-values  yTnD; 

,,  .     t   ~YLPR 
PIII(  }  =  6 


If  p   (t)  >  1,  it  is  replaced  by  1 


For  each  estimate,  the  average  bias 


.   500  . 

B  -  e^o  I   [Pi(t)  ~  P(t)]  {H'H) 

i  =  1 


and  the  relative  mean  square  error 
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Rel  MSE  = 


500 


500 

I 
1-1 


p.(t)   -   p(t) 

pltT 


T     2 


(4.5) 


are   computed  where 


p(t)    =le"(lJ 


-  Dt 


(4.6) 


The  results  appear   in  Table   3. 


TABLE   3 
Mean  Bias   and  Mean  Relative  Square  Errors 
for   Estimates   of   p(t). 


Distribution 

Exponential 

Exponential 

Exponential 

-  =.6,  T-1 
p(1)  =  .3081 

-  =.7,  T=2 
p(2)  =  .2971 

1  =.8,  T=3 
p(3)  =  -3779 

Estimate 

B    Rel  MSE 

B   Rel  MSE 

B   Rel  MSE 

(S.E)   (S.E.) 

(S.E)   (S.E.) 

(S.E)   (S.E.) 

I 

.004    .138 

.024    .328 

.005    .326 

(.005)  (.009) 

(.008)  (.026) 

(.010)  (.020) 

II 

.044    .216 

.060    .409 

.054    .436 

(.006)  (.020) 

(.008)  (.032) 

(.011)  (.031) 

III 

.0138   .137 

.040   .354 

.034   .397 

(.005)  (.009) 

(.008)  (.030) 

(.011)  (.029) 
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The  entries   for   estimate  II    in  Table   3  suggest   that  jackknif ing  <  and 
then  computing  the   estimate  of   the   probability  using  the  jackknifed  estimate 
of   k   increases   the   bias   and  relative  mean  square  error   of   the  estimate  over 
p    (t)    the  straight-forward  estimate.      The  entries  for   estimate   III   suggest 
the  jackknif ing  ln[t   -    In  p(t)]  gives   about   the   same  relative  mean  square 
error   as    the  original   estimate   pT(t)   but    increases    the   bias   somewhat. 

A   simulation  experiment   was   conducted  to  investigate  the   performance  of 
confidence   intervals  obtained   by  jackknif ing  ln[t   -   In  p(t)].      For   each 
simulation  replication,   the  average  and   variance  of  the  pseudo-values  are 
computed; 


10 


'j   "      10     *     y(j) 


1    10  -2 

°2J  -   9    \    (y(j)"y) 


An   80$  confidence   interval   for  ln[t  -   In  p(t)]    is 


[LPL.LPV]    =  Yj±   (1.383) 


10 


where    1.383   is   the   80%  point  for   a  Student's   t-distribution  with   9   degrees 
of   freedom.      The   limits   of   the   interval   are   inversely  transformed  to  give   a 
confidence   interval   for   p(t). 
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and 


di  r«.  LPVi 

PL  =  exp{t-e        } 


PV  =  exp{t  -  eLPL}. 


If  a  confidence   limit   exceeds    1,    the   limit   is  set   equal    to   1. 

Results   of   the   confidence   interval   simulation  appear   in  Table   4. 
Reported  are  the  number   of   the   500   intervals   that   cover   the   true   value 
p(t)    (C);    the  number   of   the   500   intervals  for  which  the  upper  limit  PU   is 
below  p(t)    (B);    and  the  number   of   the   500   intervals  such  that   the   lower 
limit   PL   is  above  the  true   value    (H).      The  average  width  of  the   confidence 
interval    is  also  given. 
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Table  4. 
80J  Confidence  Intervals  for  p(t 


Coverage 


Average  Width 


B 

C 

H 

(standard  error ) 

{%) 

{%) 

(%) 

Exponential 

-  =.6,    T-1 

56 

415 

29 

•  325 

p(1)    =    .3081 

(11.2) 

(83) 

(5.8) 

(.006) 

Exponential 

-  =  .7,    T=2 

52 

408 

40 

.477 

p(2)    =    .2971 

(10.4) 

(81.6) 

(8) 

(.009) 

Exponential 

-  =.8,    T=3 

54 

412 

34 

.597 

p(3)    =    -3779 

(10.8) 

(82.4) 

(6.8) 

(.010) 

The  coverage  of  the  confidence  intervals  is  within 


.80  ±  (1.96)  /^'81qq2)    =  -80  ±  -°35  =  [-765,  .835].   The  width  of  the 
interval  increases  as  —  increases. 


4.3  Simulation  Results  for  Estimating  p(t)  for  Mixed  Exponential  Service 

Times 

In  this  subsection,  results  are  given  of  a  simulation  experiment  to 
estimate  p(t)  =  P{W>t}  in  the  case  of  mixed  exponential  service  times.   For 
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each  replication  50  random  numbers  are  generated  from  a  mixed  exponential 
distribution  with 


1  "V   1  "M2t 
P{S>t}  =  ^e     +  *e  ,         t>0. 


The  estimate  ln[t   -   lnp(t)]    is  jackknifed  with  ten  subgroups   as   before  where 
p(t)    is  given   by    (3.5)   and  80%  confidence   intervals  are  constructed.      Each 
simulation  has   500  replications. 

The   asymptotic  distribution  of   the  virtual  waiting  time,  W,   can  be 
found   by  inverting  the  moment  generating  function   (1.1);    it  is  a  mixture  of 
two  exponential   random   variables. 

Table   5  reports  results   of   three  simulations. 


Case   I:  -=.30  -=.90,  T=1, 

y1  U2 


E(S)    =   0.6,       P{W>T}    =    .3^04; 


Case   II:  -     =   .35  -     =   1.05,  T  -   2, 

y1  M2 


E(S)    =  0.7,      P{W>T}    =    .3459; 


Case   III:  -     =   .40  -     =   1.20,  T  -  3, 

M1  »2 


E(S)    =  0.8,      P{W>T}    =    .4333; 
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For  each  simulation,  the  mean  bias  and  relative  mean  square  error  (4.4) 
and  (4.5)  are  computed  with  p(t)  =  P{W>t}.  Confidence  interval  coverage  and 
average  widths  are  also  computed. 

Both  cases  II  and  III  each  had  one  replication  for  which  using  all  the 
data  to  compute  the  estimate  p(t)  of  (3.5)  resulted  in  a  value  larger  than 
1 ;  these  two  replications  are  not  counted  in  the  summary  statistics  in  Table 
5. 

Comparison  of  the  mean  biases  and  mean  relative  square  errors  of  Table 
5  with  those  of  estimate  III  in  Table  3  shows  that  they  are  about  the  same 
for  both  service  time  distributions.   The  coverage  of  the  confidence 
intervals  in  Table  5  is  again  within  [.765,  .835].   However,  the  average 
length  of  the  confidence  intervals  for  the  mixed  exponential  cases  are 
larger  than  those  reported  in  Table  4  for  the  exponential  distributions. 
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Table  5 

Results  for  Estimates  of  P{W>t} 

for  an  M/G/1  Queue  with  Mixed 

Exponential  Service  Times 


Distribution 

Estimate 

Confidence 
Coverage 

Intervals 

Width 

Bias  Rel  MSE 

B     (C) 

H 

Average 

Case 

(S.E)   (S.E.) 

%           (%) 

% 

(S.E.) 

I 

.024    .172 

58   (407) 

35 

.392 

(.006)  (.013) 

11.6  (81.4) 

7 

(.007) 

II 

.040    .384 

63  (408) 

28 

.500 

(.009)  (.028) 

12.6  (81  .4) 

5.6 

(.010) 

III 

.022    .350 

54   (411) 

34 

.637 

(.011)  (.021) 

10.8  (82.4) 

6.8 

(.011) 

4.4  Simulation  Results  for  Estimating  P{W>t}  for  Gamma  Distributed  Service 

Times 

In  this  subsection,   results   are  given  of   a  simulation  experiment   to 
•estimate  p(t)   =   P{W>t}    in  the  case  of  gamma  service  times.      Each  experiment 
has   500   replications.      For   each  replication,    50   service   times   are  generated 
having  the   distribution  of   the  sum  of  two  exponential   random  variables  each 
having  mean  -.      The  estimate  ln[t   -   In  p(t)]    is  jackknifed  with  ten 
subgroups   as   before  where   p(t)    is  given   by    (3.5)   and  80%  confidence 
intervals   are   constructed. 
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The   asymptotic   distribution  of   the   virtual   waiting  time,   W,   can   be 
found  by   inverting  the  moment   generating  function    (1.1);    it   is  a  mixture  of 
two  exponential   random   variables. 

Table  6  reports   results   of   three  simulations. 


Case   I:  -  =    . 30  T  =    1  , 


E(S)    =   0.6,       P{W>T}    =    .2508; 


Case   II:  -  =    .35  T  =   2, 


E(S)    =   0.7,       P{W>T}    =    .2232; 


Case   III:  -  =    .40  T  =   3, 


E(S)    =   0.8,       P{W>T}    =    .2950; 
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Table   6 

Results   for   Estimates   of   P{W>t} 

for   an  M/G/1    Queue  with  Gamma 

Service  Times 


Distribution 

Estimate 

C< 

Confidence 
average 

Intervals 

Width 

Bias  Rel  MSE 

B 

(C)    H 

Average 

Case 

(S.E)   (S.E.) 

% 

(*)    % 

(S.E.) 

I 

..007    .110 

61 

(401)   38 

.227 

(.004)  (.007) 

12.2 

(80.2)  7.6 

(.003) 

II 

.0167   .299 

61 

(397)   42 

.328 

(.005)  (.026) 

12.2 

(79.4)  8.4 

(.007) 

III 

.0251   .390 

62 

(407)   31 

.490 

(.008)  (.034) 

12.4 

(81.4)  6.2 

(.010) 

The  mean  bias  and  relative  mean  square  error  are  about  the  same  as  for 
the  exponential   and  mixed  exponential   service   time  distributions.      The 
confidence   interval    coverage  is  once  again  within    [.765.    .835].      The  average 
width  of   the   confidence   intervals   is   smaller   than  the  widths   in  the 
exponential  and  mixed  exponential   cases.     This  is  to  be  expected  since  the 
gamma  distribution  has   a  shorter   tail    than   the  exponential   or  mixed 
exponential   distribution. 
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